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1 Abstract 


The first steps in the neural processing of sound are located in the auditory nerve 
and in the cochlear nuclei. To model the signal processing efficiently, we propose a 
simple mathematical tool that takes the minute timing of the system into account. 
In contrast to the situation in the cortex, the number of connections between neurons 
in auditory periphery is comparatively low. This gives way to an accurate modeling 
of the connectivity of the neuronal network. The timing is the all important feature 
in the peripheral neuronal auditory pathway. The primary auditory neurons e.g. 
phase lock to periodic sounds with important interactions with respect to both the 
refractory periods of the neurons and to the time delays caused by traveling times 
along the basilar membrane or through a synaptic connection. 

The mathematical tools provide a solid basis to build models for peripheral 
auditory processes. In particular, we study carefully a large class of refractory 
neurons, hud analytical formulas for the spiking activity, and prove that refractory 
neurons respond to periodic signals by asymptotically periodic output. The methods 
rely on the theory of positive operators and give a numerical scheme for finding 
fixed points to an integral operator with geometric convergence rate. In addition, 
we consider a perfect integrator neuron, mathematically equivalent to randomized 
random walk, where the random walk is bounded from below, and solve the first 
passage time problem using continuous time Markov chain techniques. Our method 
leads to ordinary differential equations that are linear. The dynamical behavior can 
thus be described by classical methods. 

In an accompanying paper we set up the simulation framework as a counterpart 
to the present mathematical model. By suitably adjusting the few parameters in 
the model it is possible to reproduce the basic patterns of neural activity. 


2 Introduction 

In this paper we try to develop a mathematical tool that allows to describe signal 
analysis in the auditory pathway of the brain. The incoming signals from the in¬ 
ner ear are processed in parallel pathways up to the inferior colliculus. Essentially 
starting from this neuronal center information from other parts of the brain are com¬ 
bined with the auditory components. A time dependent pattern of neural activity 
then emerges in the auditory cortex. Signal processing along the auditory pathway 
is complex. There is a wealth of detailed experimental information available and 
many of the essential features have been described precisely [19]. To a large ex¬ 
tent the activity of specihc neurons has been classihed and some information on the 
topological ordering of different groups of neurons is available. The challenge is to 
set up a global picture and to understand how the different components combine on 
the signal processing level. 

The approach taken here is to describe the neuronal activity with densities that 
have a probabilistic interpretation. The neuronal nuclei are then pictured as trans¬ 
forming these densities in a specific way, adapted to the physiological significance of 
these nuclei. Densities give a less precise description of neurons than models based 
on dynamical systems would admit. Yet densities should provide a more efficient 
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tool in the treatment of the interaction of the neuronal nuclei and they should facil¬ 
itate keeping track of the patterns of neural activity. They are specifically designed 
towards time and intensity coding. 

We land at the area of firing rate models (see chapter 11 of [1]) and choose a 
Poisson type model that allows to represent the firing properties of neurons on the 
basis of densities. We hrst analyze the signal processing of a single neuron with 
respect to the refractory properties. Such effects have already be studied previously 

in [3]. 

Our method is based on stochastic results from renewal theory. With it, any 
kind of refractory period can be handled. 

For periodic signals we show that any refractory neuron responds to a periodic 
input by an asymptotically periodic output. 

To model the combined activity of a network of refractory neurons, we use perfect 
integrate-and-fire with stochastic input and bounded paths. It turns out that a 
continuous time hnite state Markov chain model (CTMC) is suitable to the analysis. 
In particular, the synaptic and transmission delays are easily incorporated. 

This approach circumvents Stein’s model and stochastic differential equations. 
Yet still, it seems to capture the essential features observed with the neurons in 
the auditory pathway. Our model leads to an ordinary differential equation that is 
linear. Explicit solutions can then be obtained through standard techniques in linear 
algebra and differential equations. We e.g. get an easy derivation for the classical 
formulas for perfect integrate-and-hre neurons receiving only excitatory inputs. 

Our mathematical approach is amenable to modeling and relates to simulations. 
In a companion paper this will be investigated in detail. 

In order to focus our work, we have asked for transparency and simplicity of the 
model and our efforts are directed towards using as few parameters as possible. We 
take into account the hne timing parameters like the refractory properties of the 
neuron and the synaptic and transmission delays. We also include the spontaneous 
activity of the neuron and the local architecture of the network. We omit any neu¬ 
rochemical variables and hence the model is at best only phenomenological. 


3 Stochastic processes for auditory periphery 

A neuron receives spike trains through its dendritic tree and emits a spike whenever 
the membrane potential of the cell reaches a cell specihc threshold value. The spikes 
are narrowly supported in time and have uniform shapes. Thus the mathematical 
theory of point processes is widely used for modeling the neurons statistical behavior. 
In particular the problems related to the neuronal ensembles can be reformulated in 
terms of queueing theory. 

There are models for single neuron behavior of varying detail. The most detailed 
ones are based on the physical model by Hodgkin and Huxley and contain several 
parameters related to the chemical and physiological properties of the neuron. From 
the mathematical point of view, this model describes the voltage of the membrane 
in terms of a dynamical system, i.e. a system of ordinary differential equations. The 
simplest model is the integrate-and-hre (IF) and a more sophisticated is the leaky 
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IF (LIF), which breaks the problem of neuronal firing in two parts: the time before 
the neuron fires is modeled with differential equations and once the solution reaches 
critical level, the neuron fires and then rests before the process starts anew. 

The classical models are deterministic. Stochastic nature of the neural activity 
can be incorporated in several different ways. The Stein model [15] is a stochastic 
differential equation model for the sub-threshold membrane voltage evolution. It 
assumes excitatory and inhibitory Poisson inputs and an exponential decay of the 
membrane potential. This is still a very detailed model and correspondingly, explicit 
calculations are difficult, e.g. an explicit formula for the first passage time of the 
system with constant intensity Poisson inputs is not known (see Sacerdote, Giraudo 
2011). If the leakage term is neglected, the resulting model is randomized random 
walk (RRW) for which the first passage time can be derived - assuming constant 
input - via Laplace transforms [18]. However, there is no formulas available for 
non-constant stimuli. 

In the present approach we take every incoming signal to be inherently stochas¬ 
tic by assuming always an explicit stochastic process, which comes with a certain 
intensity depending on the (acoustic) stimulus and previous processing steps only. 
The mathematical expression for this intensity is a non-negative time dependent 
density s{t). We assume, that the neurons activity is completely described by the 
nature of the stochastic process (e.g. Poisson, Gamma type neuron), by the density 
s{t) that itself is derived from the input signals from other neurons and by the time 
of the previous spike of the neuron itself. To easily combine the activities of the 
different neurons, we only compute the output intensity of the spiking and assume 
that when several of these intensities are combined in a next neuron, the pooled in¬ 
coming spikes look like a Poisson process. In the case of constant intensity processes 
which are sufficiently regular, the sum of the processes approaches Poisson process 
(see [2]: Proposition 11.2.VI). 

We are interested in two types of abstract neurons: primary neurons are directly 
stimulated by a continuous variable like the concentration of neurotransmitter in 
the inner hair cell auditory nerve fiber complex - or a more abstract variable like a 
probability density; integrating neurons take electric spike trains from other neurons 
as input and process the spike trains according to a rule where the excitatory input 
brings the neuron closer to firing and inhibitory spikes push the neuron away from 
emitting a spike. The rules are motivated by the physiological properties of the cells. 
We try to characterize the quantitative behavior of these abstract neurons based on 
first principles. 


4 The refractory neuron 

4.1 The definition of the refractory neuron 

An abstract refractory neuron, ur, is taken to be a stochastic process determined 
by an inhomogeneous, non-negative function s{t) and a homogeneous function r{t). 
The inhomogeneous function represents the external stimulation of the neuron and 
the homogeneous function describes the internal dynamics of the neuron. As an 
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example, s{t) could describe the amount of the neurotransmitter in the inner hair 
cell, and r{t) the refractory properties of the auditory nerve hber. More precisely 
we dehne the probability of the neuron to emit a spike inhnitesimally through 

P{ti > t + h\ti > t) = 1 — s{t)r{t — to)h + 0{h‘^) (1) 


where to is the moment at which the neuron emitted the previous spike before ti 
and is an error term with 0{h)/h —)■ 0 as h —?• 0. We assume throughout 

that s{t) is a non-negative locally integrable function, r{t) = 0 for all negative t, 
and r{t) < r{u) for all t < u. Furthermore we assume that the integral 

s(t)r(t — x)dt 

is inhnite for any x > 0. 

To determine the probability density of the bring we bx f and divide the 
interval [to, t] into n equal parts with nh = t — to, aj = to + (j — l)h and bj = aj + h, 
j = 1,... ,n, to have 

n 

[to,t] = [j[aj,bj]. 
i=i 

By the debnition of the conditional probability we have 



P{ti > t\ti > to) 


Pjh > t) 

P{ti > to) 

Pjh > bn) Pjti > bi) 

P{ti > an) ” ' P{ti > ai) 

n 

> bj\ti > aj). 


Taking logarithm, using the assumption © and linearizing the right hand side 
logarithms gives 

n 

logP(fi > t\ti > to) = lim y —hs{aj)r{aj — to) 

n —>’00 ' ^ 

which can be interpreted as an integral to give 


P{t\ > t\ti > to) = exp 


s{u)r{u — to)du 


'to 


Let us debne now the transition probability p{x, t) by diberentiating the conditional 
probability above with respect to t and setting to = x to obtain 


p{x, t) = s{t)r{t — x) exp 


{u)r{u 



( 2 ) 


Then p{x, t) is the probability density that the brst bring after bring aXto = x occurs 
at t. We set 

p{x,t) = 0 for t < X. (3) 
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Clearljfl 


p{x, t)dt = 1 


(4) 


If the “first” firing time has a distribution po, then the next brings are given recur¬ 
sively by 

Pk{t)= pk-i{x)p{x,t)dx, 

J —OO 

for fc = 1, 2,.... The consecutive bring times form a point process 
Analogously, the transition probabilities for the fcth bring are 

Pk{x,t)= pk-i{x,z)p{z,t)dz, 

J —OO 

Similar derivations can be found in [8]: 5.2.3, and m- pp. 4-5. 


4.2 The output of the refractory neuron 

Given either a refractory neuron or a integrate-and-bre neuron as a model for the 
“law” of the neuron, we debne the output rate of the neuron. Let T = ... } 

be a point process. The attached counting process is debned by 

N{t) = {i> 0\U <t< U+i}, 

and we call M{t) = E{N{t)) the expected number of neural spikes emitted up to 
moment t. If M is diberentiable, we debne 

_ dM{t) 
dt 

to be the instantaneous bring rate. This is the time varying output rate of the 
neuron’s activity. In the physiological measurements, the peri stimulus time his¬ 
togram (PSTH) corresponds to the instantaneous bring rate. On the other hand, in 
stochastic analysis M{t) is called the renewal function. 

Mathematically, in the simplibed model treated in this paper, a simple neuron is 
the transformation of incoming rate functions s_ and s+ to an output rate function 
lit). 


4.3 Instantaneous firing rate 

Given a locally integrable non-negative function s{t) and an initial probability dis¬ 
tribution pidt) = poit)dt for the brst spike, we compute the instantaneous bring 
rate /(f). 

The instantaneous bring rate I{t) of a neuron at time t corresponds to the joint 
probability density function pix, t) of the neuron to bre at moment t given it bred 
at X and the initial probability measure p due to the history of the system up to 

^We often use the shorthand J for the definite integral and likewise we abbreviate = 
, if not otherwise stated or clear from the context. 

— OO ■ 
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a starting moment t = 0 (stimulus onset). We have (justified by [9]: Chapter 7.7: 
Theorem 4) 

/*t2 pt poo 

m = E/ ■■■/ / /i(s)p(s, ti)... p{tk, t)dti... dtkds (5) 

ti=0 ik—tk — 1 *^0 

oo 

= '^Pk{t) (6) 

k=l 

The definition of instantaneous firing rate is related to the repeated measurements 
data gathering procedure, where the activity of single neurons is recorded while 
the auditory system receives acoustic input. It is tacitly assumed that between the 
stimuli there is enough time for the neural system to get back to a equilibrium state 
where only the spontaneous activity of neurons persists. If the spontaneous activity 
is assumed to be a homogeneous Poisson process, then the corresponding measure 
yU would be just the exponential function 

fi{s) = Xe~^^ds. 

In a mathematically simpler situation, we can assume that the neuron fired at f = 0, 
at the onset of the stimulus. This corresponds to the choice 

/i(s) = 5(s), 

where 6 is the Dirac measure concentrated at the origin. 

The case with no refractory period ( r{t) = Xoi^) ) corresponds to the classical 
Poisson model. In this situation, the instantaneous firing rate gives back the original 
intensity 


Pm{x;t) = dtXx{t)s{t)exp{- s{u)du) 


"Vs 


''y2 


s{ym-i)dy^-i... / s{y2)dy2 / s{yi)dyi 


Upon setting w{t) = J^s{u)du, s{t) = one obtains 


i-ys 


ry2 


s{y,n-i)dym-i-- / s{y2)dy2 / s{yi)dyi = 


{m — 1)! 


and the transition probabilities have the form 

Pn,{x]t) = Xxit) s{t)e-^d) - (7) 

(m — 1)! 

The instantaneous firing rate of a neuron at time t, given that the neuron fired 
at time x, is 


= Xx{t)s{t)e = Xx{t)s{t). (8) 

m=l 
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We thus obtain the original intensity s and observe that the instantaneous bring 
rate does not depend on the past. Hence, 


lim g(x;t) = s(t) 


for all t and the convergence is uniform on any bnite interval. In the following 
sections neurons with non-trivial refractory period will be considered. 

4.4 Densities 

At the level of single neurons, we assume that the activity of a neuron is determined 
by the incoming densities Si{t) > 0 and a set of parameters: spontaneous activity of 
the neuron a, connection strength a;*, connection delay Tj, refractory function r and 
bring threshold 'd. We brst divide the analysis of the neuron into two independent 
steps: input analysis and output generation. The brst results in a single pooled 
activity, a Poisson process that we approximate by the integrated input density s. 
The output generation depends on s and r, the refractory component, as well as on 
a and -d. The spontaneous activity a can be subsumed in the calculation of s{t). 

The inbnitesimal probability density that determines the instantaneous bring 
rate of the output is a product of the integrated input density s(t), and the refractory 
function r{t — x), translated to the position x of the previous bring time (see |12] . 
[7], and [13], the details are given in the next section). 

A typical example for r is given by 



0 if M < Pa, 
1 - exp(-u/pR) ifa > Pa, 


(9) 


where the constants pA and pR are absolute and relative refractory periods respec¬ 
tively. We assume that r is always monotone increasing, grows at most polynomially, 
and r{u) = 0 for all u < 0. We do not assume continuity for r. 

For the incoming densities we consider two cases: either we have an auditory nerve 
bber, for which a simple inner hair cell model will then provide the input density, 
or the input comes from a collection of neurons with specibc outputs. In the latter 
case, we represent the activities of the neurons that connect over the dendritic tree 
by densities Si,i G I. 

The (excitatory) inputs are pooled by taking the sum of the densities s*, weighted 
with the strength uji of the connections, and taking into account the delays of 
each path. 



( 10 ) 


where I is the set of indices for the (excitatory) connections. 

For the refractory neuron only excitatory inputs are taken into account. The 
situation of mixed, excitatory and inhibitory inputs, will be studied in the context 
of the integrate-and-bre neuron. 


The pooling of the densities is motivated by the following addition property 
Assume that the neuron Y fires whenever any of the two input neurons Xi or X 2 fire. 
If the firing densities for Xi and X 2 are si and S 2 respectively, then the resulting 
output density of the neuron Y is s = Si + S 2 
This statement fails, if refractory periods are involved. 

The proof is a direct consequence of the exponential law: 


P(Ti^ > t\To <x) = 1 - P{Tf^ < < x)P{Tf^ < < x) 




Xi 


- 1 X 2 


-1X2 


= 1 — exp I — / si{u)du exp I — / S2{u)du 


= 1 — exp 


( si ( m ) + S2{u))du ) . 


Hence, Y can be obtained by summing the densities of Xi and X 2 . 


4.5 Periodic signals 

The transition probability (see equation ([2])) associated to the neuron and to the 
density s{t) is 

p{x,t) = r{t - 

For the periodic case it is assumed that s is 1-periodic. As a consequence 

p{x + m,t + m) = p{x, t) 


for any integer ru¬ 
in the periodic case, the instantaneous bring rate is well debned and asymptoti¬ 
cally does not depend on the initial measure /r. This is a result of Thorisson (El. 
Theorem 2). In the following we use p. = 6x the dirac measure at x. With this, the 
instantaneous bring rate of a neuron at t is given by 

00 

(l{x,t) = ^Pm{x,t) ( 11 ) 

m=l 

Under the additional condition that the factor m debned in equation flT^ below 
is bnite, the limit 


q{t) := lim q{x, t) 

x^—oo 

exists (El. Theorem 6). 

In the discussion above, when no refractory period was involved, it was shown that 
q{t) was equal to the initial density s{t). In the present situation this will be diberent. 
Note that the function q is periodic: 

CXD 00 

q{x, f -f 1) = ^ Pm{x, t + 1 ) = Pm{x -l,t) = q{x - 1, t) 

m=l m=l 
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q{t+l)= lim q{x — l,t) = q{t) 

X—>• —oo 

The operator T on locally integrable periodic functions is defined as 

Tf{t) = j p{x,t)f{x)dx 

and its dual on by 

T*g{x) = J p{x,t)g{t)dt 
The fact is then that q is invariant under T 

= / p{x,t)q{x)dx = / p{x,t) lim q{z,x)dx 

I I z^—oo 


/ oo 

p{x, t) y 

m=l 

CXD 

= lim 

Z^ — OO ^ ^ 

m=l 

= Mm {q{z,t)-p{z,t)) 


Pm[z,x) 


= g{t) 


The operator T can be periodized. 


Tf{t) = j p{x,t)f{x)dx 

= / p{x + m,t)f{x + m)dx 

= f p{x + m, t)f{x)dx 


Set 


m] 


p{x,t) = E p{x + m, t) = '^^p{x, t + 

m m 

then the operator T on periodic functions reduces to 

Tfif) ■= [ p{x,t)f{x)dx 
Jo 

on the space T^[0,1]. Its kernel is periodic in both variables and satisfies 
/ p{x,t)dt = / '^^p{x + m,t)dt = / '^^p{x,t + m)dt = / p{x,t)dt 

JO JO ^ ^ J 

Furthermore, q considered as a function in T^[0,1] is invariant 

git) = Tqit) = fqit) 
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In the case of a 1-periodic infinitesimal density function it is possible to start 
with the probability distribution Q{A) that the neuron has hred at to ^ ^ mod 1. 
The set A is a subset of [0,1) and Q([0,1)) = 1. Assume then that the neuron has 
hred at time tg = x within the period. 

The initial probability Q is transformed into the probability 



p{x, t)dQ{x) 


If Q is mapped onto itself under this transformation, then it can be represented 
by a density qdt 

Q{A) = [ q{t)dt 
Ja 

that satishes the hxed point equation 


q{t) = / p{x,t)q{x)dx. 


The output density q{t) of the neuron driven by the periodic input density s{t) is a 
hxed point of T and hence of the periodized operator T. In general it will not be 
normalized whereas the hxed point q above is normalized by q{t)dt = 1. 

Theorem 4.1. (Thorisson) The integral operator 


Tq{t)= / p{x,t)q{x)dx. 


has a unique normalized fixed point q G L^([0,1]). Furthermore for every non¬ 
negative qo G T^([0,1]) with ||go|| = 1 sequence qn = Tqn-i,n = 1,2,... converges 
exponentially to q in T^([0,1]). 


This result is contained in [HI- An independent proof will be given in Appendix 
2 fCorollary 16.4p . 

The uniqueness statement tells us in particular that the outgoing inhnitesimal 
probability distribution q of the neuron - provided it exists - is a multiple of the 
normalized distribution q. 

The expectation time for ti to occur on the condition that to = x is 


jtpOXdt 


The expected delay is thus 


E{x) 




j(t-.)pO,t)dt 


Theorem 4.2. (Thorisson) Assume that 


m = 


E{x)q{x)dx < oo 


( 12 ) 
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where q is the normalized fixed point above. Then the outgoing infinitesimal proba¬ 
bility distribution approaches the periodic density 

1 „ 


r-n+l 


This entails that lim 


|g(a;, t) — q{t)\dt = 0 for every x. 


This is proved in HZI, Theorem 3. As an example take the case of a nenron 
with an absolnte refractory period g and a constant input density s{t) = A. Here, 
q{t) = 1 is a fixed point for T. The expected average delay is 


•*1 noo 


E{x)dx = dx {t — x)Ae ^^~^P">^dt 
Jo J x-\-p 


-1 


The outgoing density is thus 


= p + A 


QiO = 


p +A-i 

in accordance with the result in Appendix 1. 


4.6 Constant stimuli with general refractory structure 

We consider an important special case which corresponds to a neuron with sponta¬ 
neous firing and a short refractory time. Typically the spontaneous firing rate is at 
most 100 Hz while the refractory time is at least 0.7 ms. We model the situation by 
taking a constant density. 

In the special case s(t) = A, A > 0, the transition probability 

is a function of t — x that will be written as p{t — x). 

The transition probabilities for the k-th bring are then obtained by convolution 

Pk{t-x) = Jpk-i{t - z)p{z - x)dz 

= pk-i*p{t - x) = p*^{t - x) 


Provided the absolute refractory period is positive, the instantaneous bring rate 

OO 

q{t -x) = '^Pmit - x) 

m=l 


with system start at x has only bnitely many terms. The instantaneous bring rate 
of the neuron is then given by 


q 


lim y^Vkit-x) 

-A^ — OO 

1 

OD 

lim 
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In this situation one actually need not invoke Thorisson’s theory, as the classical 
Blackwell renewal Theorem [5] verihes that the previous limit exists and is constant. 
In the case of a pure absolute refractory period p > 0 we already noted that 

A 


Another proof that uses just simple Fourier analysis in included in Appendix 1. 

5 The integrate—and—fire neuron: a simple dis¬ 
crete model 

In this section we derive the instantaneous bring rate and the interspike interval 
distributions for the integrate-and-bre neuron. We analyze the example of constant 
inputs and bnally show that for every periodic stimulus there exists a unique periodic 
equilibrium density which the neuron’s instantaneous bring rate approaches. This 
simple model does not allow the inclusion of a refractory period. 

5.1 Integrate-and-fire neuron 

An integrate-and-bre neuron (IF neuron), riip is debned through inhibitory and 
excitatory incoming spike trains. The model described here is perfect in the sense 
that the ebect of the incoming spikes does not decay over time. Moreover, the model 
has bounded paths since the neuron is not allowed to have inbnitely big membrane 
potential values. This is done by modeling the membrane potential in an abstract 
way as states in a bnite system, where each state characterizes how many excitatory 
incoming spikes are needed at least before the neuron can emit a spike. 

More precisely, let S-(t) and s+(t) be the intensities of two independent, inho¬ 
mogeneous Poisson processes. The IF neuron has K > 1 possible states 1,... ,K. 
When a spike arrives from the inhibitory process, the IF neuron moves from state i 
to state i + 1 unless already at the lowest state K. Similarly, when a spike arrives 
from the excitatory process, the neuron moves up from the state itoi — 1, ifi>l 
and moves from state 1 to K otherwise. This transition is called resetting and the 
IF neuron emits a spike during the transition. These transition times form a point 
process Tip = {to, ti, • • • }• 

5.2 Densities 

We assume that the activity of the neuron is determined by the incoming densities 
Si{t) > 0 and a set of parameters: spontaneous activity of the neuron a, connection 
strength cuj,connection delay Tj and bring threshold d. We divide the analysis of the 
neuron into two independent steps: input analysis and output generation. The brst 
results in a single pooled activity, a Poisson process that we approximate by the 
input density s. The output generation depends on s as well as on a and 'd. The 
spontaneous activity a can be subsumed in the calculation of s{t). 

The incoming densities s*, i & I. represent the activities of the neurons that connect 
over the dendritic tree. The excitatory and inhibitory inputs are pooled separately 
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by taking the sum of the densities Sj, weighted with the strength uji of the connec¬ 
tions, and taking into account the delays Xj of each path, 


(t) = '^UiSi{t - Ti), 
iei+ 

(13) 

(t) = - Ti), 

iei- 

(14) 


where G I and I- G I are the subsets of indices for excitatory and inhibitory 
connections respectively. The probability of an incoming spike is modeled by the 
density s{t) = + s~{t). 


5.3 The infinitesimal generator 

A probability vector v = {vi,V 2 , ■■■,Vn) is a vector with non-negative components 
that add up to 1: n* > 0 for i = and = 0. A linear mapping that 

maps probability vectors into probability vectors is called a probability mapping. 
Expressed by a matrix A = (aij) in standard coordinates it is characterized by the 
condition that its column vectors are all probability vectors. If the linear differential 
equation 

v' = Qv 

generates a flow of probability mappings, then Q is called an inhnitesimal generator 
for probability mappings. It is well-known that the inhnitesimal generators Q can 
be characterized by the property that they can be written in the form Q = {A — I)s, 
with A a probability matrix and s > 0 is a scalar. In other words, the diagonal- 
elements of Q are non-positive, other elements non-negative, and each column sums 
to zero. Note that the how of the diherential equation maps the positive cone 
= {v : Vi > 0,i = 1, ...n} into itself. In fact, on the boundary of the cone, the 
vector held Qv points into the cone or is tangent to the cone. 

The inhnitesimal generator of the process described in subsection 15.21 is given by 
the matrix Q{t) = {A{t) — I)s{t) with a matrix A{t) of the form 


0 p{t) 



q{t) 0 

pit) 



(l{t) 

0 pit) 

pit) 


qit) qit) 


with p{t) = ^ and q{t) = 

The time development of the system is fully determined by the matrix diherential 
equation 

v'{t) = Q{t)v{t), 

where n is a time-dependent probability vector. 
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A first order linear matrix equation has always a unique solution with prescribed 
initial condition v(x) = w. It can be expressed by the Peano-Baker series in terms 
of iterated integrals, see El. i. ra. 


v{t) = x)w 


where 


^{X,t) = I + / Q{Ti)dTi + 


t fTl 



Q{Ti)Q{T2)dT2dTi + ... 


5.4 The periodic case 

Theorem 5.1. Assume that A{t) is a 1-periodic probability matrix that is piecewise 
continuous in t, and that s{t) a 1-periodic function. Then the equation 

v' = {A{t) — I)s{t)v (15) 

has a 1-periodic solution. 


This is a consequence of Floquet theory: 

For all t, the adjoint matrix A{t)* has eigenvector (l,l,...l) with eigenvalue 1. 
Therefore the adjoint equation 

-y-= {A{t) - iys{t)y (16) 

has (1,1, ...1) as a constant and hence periodic solution. Following [TT] p. 94, the 
space L of periodic solutions of (1) has the same dimension as the space L* of 
periodic solutions of (2). Therefore the space L has dimension at least one. 

Theorem 5.2. In this situation, the (normalized) periodic solution is unique. 

The proof is given in Appendix 2. 

\i V = {vi,V 2 , ...,Vn) is the solution of the equation, then the output density of 
the neuron is given by the density 

p{t)vi{t)s{t) = s+{t)vi{t) 

that controls the passage from the hrst to the n-th level of the neuron. 


5.5 Interspike interval distribution 


Given a vector w representing the distribution of probability over the states at the 
initial time x, we determine the transition probability p given the time-dependent 
Poisson processes N+ and iV_ as input. By adding an absorbing state in the diagram 
77 we hnd the inhnitesimal generator matrix 


Qoit) 


0 s+{t) 

0 —s{t) s+(t) 

0 S-{t) —s{t) s+{t) 

0 S-{t) “^(t) s+(f) 

0 S-{t) —s+{t) 
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where s(t) = s+(t) + S-(t). The time development of the finite state system up to 
the next firing is the solution to the differential equation 


v'it) = Qoit)v{t) 


(17) 


The solution v(t) = {vo{t),vi(t), ...Vn(t)) of the differential equation with initial 
condition v{x) = e„ is the probability distribution over the states, given that the 
neuron bred at time to = x (at this time the neuron is at state n). The component 
Vo{t) is the probability that the next firing occurred in the interval {x,t]. The 
transition probability is therefore given by 



( 18 ) 


The output density s* (that depends on x and on t) of the neuron can then be 
calculated as the instantaneous bring rate 



1 


(cf. formula (10)). 


5.6 Refractory period for the integrate and fire neuron 

The model for the integrate and bre neuron can also be considered as a renewal 
process with time dependent input. The process depends on the last bring time x 
and on the excitatory and inhibitory inputs densities s+ and s_. As in section 4.4 
it is then possible to incorporate a refractory period in the model. The refractory 
function r{t — x), translated to the position x has to be multiplied with the input 
densities. The diberential equation v'{t) = Q{t)v{t) is then solved with s+(t) re¬ 
placed by s+(t)r(t — x), S-(t) by S-(t)r(t — x) and consequently s(t) by s(t)r(t — x). 
The transition probability is again given by (1T8|1 . with vi{t) the component of the 
solution vector of the modibed diberential equation. 

5.7 Constant stimuli 

In general, the explicit solution of the diberential equation is difficult to bnd. This 
can already be seen by considering the simple case, where the system is time- 
independent, i.e. the excitation and the inhibition are homogeneous Poisson pro¬ 
cesses. Then the Peano-Baker series simplibes to matrix exponentiation. However, 
even this cannot be calculated explicitly for large matrices. 

Two state system 

The simplest case is the two-state system which bres at the moment t = 0. There, 
the inter spike interval density can be calculated (using Mathematica) directly and 


we have 



,2 
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Similar calculations reveal that the instantaneous firing rate with initial state 
w = (0, 0,..., 0,1) at moment t = 0 gives 

-- (1 

For larger matrices the formulas are similar but become enormous. The important 
phenomenon is that the instantaneous bring rate always approaches at an exponen¬ 
tial rate a unique equilibrium. 

Many state system with excitation only 

If inhibition is absent, then the density of the brst passage time of the integrate- 
and-bre neuron is 




X / N N exp(—s+t)(s+t)^ 

f{t) = s+{exp{Qot)w)i = - — -, 

with the choice w = (0, 0,..., 0,1). We bnd thus the well known formula for perfect 
integrate-and-bre neurons with excitation only. 

Equilibrium in the presence of inhibition 

In the equilibrium, the iF-state system satisbes Qv = 0. This leads to 

p,., Ehi(A' - 

For >> s_ the excitatory term dominates and the expectation is roughly K/s^; 
for = s_ the expectation is 


m 


K{l + K)_ 


2sj 




for << s_ the expectation blows up. 
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6 Proofs for the theorems 


Appendix 1 


Theorem 6.1. Let X be a Poisson process with a refractory time p > 0 and constant 
density s(t) = A > 0. Then the instantaneous firing rate is well defined and 


Qit) 


A 

l + Ap 
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The Fourier transform of 


P(^) = ^X{t>p}e 


is 

poo poo 

P{u) = A e-^^^e-A^-p'^dt = Ae-^^P / 

Jp Jo 

Ae-i^p 
A + ioj 

By the convolution theorem 




p*^{ijj) = A^ 


^—ikujp 

{A + iuY 


Since Pk{.x) > 0 has support on [kp, cx)) and f pk = 1, we deduce that 
|g(t)|(l + \t\)~‘^dt < oo, where q is the locally integrable function 


k=l 

and its Fourier transform satishes 

CXD 

q{u) := '^^pk{uj) convergence in S'. 

k=l 


We may perform the summation in terms of geometric series for u 0 (with locally 
uniform convergence outside the origin), and obtain 

Q-iuip 

q(u)) = - ^^—, a; 7 ^ 0. 

' ]__!_*“ _ Q—lUJp^ ' 

The function a; i—)■ -—is not integrable at the origin and hence the classical 

1+^—e P 

dehnition of the Fourier transform does not apply. However, it certainly dehnes a 
distribution in the principal value sense and below when dealing with this function 
as a distribution we assume that it is dehned as a principal value sense. Since the 
Fourier transform is unique, q will then have the form 

^ Q-iu)p 

q = d + \p.v.} - ^^—, 

^ ^ — 0—itJp ’ 

where d is a distribution supported at origin. Our aim is to show that d is a constant. 
Consider the behaviour of g — d. When uj is close to zero we have 


q{uj) — d{u) = 


A 


+ 0 ( 1 ). 


1 + Ap iu 

Since the term the numerator corresponds to a shift we neglect it for 

the moment. For big values of u 

1 A 


1 + f 


ioj + 0(1) 
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Let us write 


1 


1 + ^ - 


= R{uj) + S{uj) 


with 


R{u:) = - 


^ 4(1 iuj') 


iuj{l + Ap + ioj) 


For small values of u: we see that 


^(a;) = 

and for large values of uj we have 


A 


^ p^-2p 


S{u) = 


l + Ap 

A^{1 - p - e-^^p) 


+ 0(a;) 


+ 0{uj) 


Hence S is integrable and by the Riemann-Lebesgue lemma the inverse Fourier 
transform is uniformly continuous and vanishes at inhnity. Since 


R{u:) = 


A 


1 A^p 

+ 


1 + Ap ioj 1 + Ap 1 + Ap + ioo 


we can compute the inverse Fourier transform (recall that we need to interpret the 
distributions at the origin in the principal value sense) as 




l + Ap 


2(1 + Ap) 


sgn(t). 


But we know that q{t) is zero for t < p. Also, every distribution d supported at the 
origin has an inverse Fourier transform, which is a (complex) polynomial. Since 

q(t) = d(t) + R(t - p) + S(t - p) 

we conclude that d must be constant. Hence 

A 


d = 


2(1 + Ap) 


and 


1 . / ^ A 

q = hm q{t) = 

t—)■—oo 1 -j- j\p 


This proves the theorem. For the Fourier analysis results we refer to Hi. 


Appendix 2 

In this appendix we give a full proof of theorem 15.21 The proof is a fixed point 
argument. We show that the statement of the theorem is equivalent to hnding a 
unique hxed point of an integral operator. This operator then is shown to be positive 
and by examining the dual space, we conclude that the hxed point is attractive in 
the sense that any initial value will converge to the same hxed point geometrically 
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fast. Here we recall some classical facts on uniformly ergodic Markov chains. Let 
{S,B) be a measurable space and denote by A4{S) the Banach space of all hnite 
measures on {S,B), equipped with the total variation norm, denoted by || • ||i. 
The subset of probability measures is denoted by V{S). Let us call any linear and 
bounded operator T : A4{S) —?• A4{S) that maps positive measures to positive 
measures in an isometric way a transition operator. Especially, transition operators 
map probability measures to probability measures. The Dobrushin coefficient of 
ergodicity is dehned by the formula 


6{T) := sup 


\\Tpi - Tfi2\\i 

ii/'i-wiii 


(19) 


where the supremum is taken over all distinct ^ It is straightforward 

to check from the dehnition that for products of transition operators one has 


5{STU) < 5{T). (20) 

Theorem 6.2. Assume that the transition operator T : AA{S) AA{S) satisfies 
h{T) < 1. Then there is a unique equilibrium distribution vr G V{S) with Ttt = tt. 
Moreover, one has the exponential convergence 


TT — T’^/i||i < 2(5(T))”, for all n>l and pEV{S). 


Proof. The existence and uniqueness of tt follows by Banach’s hxed point theorem 
since the assumption 6{T) < 1 states that T is a strict contraction on the closed 
subset V{S) C M(S). Also, we may estimate for n > 1 

IItt - T>||i = llTvr - T(T"-i)p||i < 5(T)||vr - T-'VlIi, 


and the second statement follows by iteration. 


□ 


The refractory neuron 

There is a particular class of transition operators that can be expressed as integral 
operators and for which the above Lemma applies directly. In particular, this is the 
case for the operator in Theorem 14.11 in our setup. Namely, assume that the space 
(S, B) admits a reference probability measure m E P(S) and the transition operator 
T : A4(S) —>■ A4(S) acts via 

Tp,{A) = / k{x,y)m{dx)pL{dy) for any yEV{S), 

Jaxs 

where the integral kernel k : S x S ^ [0, oo) is measurable and non-negative. 
Theorem 6.3. Assume that the kernel k of the transition operator T satisfies 

k{x,y) > g{x) for all x,y E S, 

where g{x) > 0. Then 

6{T) < 1 — / g{x)m{dx). 

Js 

Especially, if the function g is not zero a.e., then 6{T) < 1 and the conclusions of 
Lemma Ih.gl apply. 
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Proof. Assume a := g(x)m(dx) > 0. Then [/ : M(S) M(S) is a transition 
operator, where 1/ has the kernel k{x,y) := a~^g{x) for x,y G S. If we dehne 
A G V{S) by setting A(A) = a~^ J^g{x)m{dx), it follows that 

Ufi = \ for any yGP{S). (21) 

The assumption on k implies that the operator T — all is positive and by applying 
it on probability measures we hence see that it satishes the norm bound 

II (T - aU) : M{S) >1(^)|| = 1 - a. 

Given any two probability measures /ii,/i 2 on S we apply fl^ to compute 

||T/ii -Tpslli = ||(T- aG)(/ii -/i2)||i < (1 - a)||hi -h2||i, 


and this proves the Theorem. 


□ 


The hrst part of the following corollary yields an alternative proof for Theorem 
14.11 For that end, given a Banach space E and Xq G i?, Xg G E' the one dimensional 
operator y i-G- (xg, y)xQ on E is denoted as usual by Xg ® Xq. 


Corollary 6.4. There is a E (0,1) such that in the situation of Theorem 4J_ it holds 
that 


11^^/- c^Uho,!) < (1 - a)*" where c := f (22) 

Jo 

for all k > 1 and all functions f G T^(0,1). 

In addition, the fixed points of the dual operator T* consist of constant functions, 
and we have 


||((T*)^ 5 ( - c||z,-(o,i) < (1 - a)'" where c := [ fq (23) 

Jo 

for all k > 1 and all functions g G L°°{0, 1). 

Proof Theorem 16.31 shows that ||T/ — ^|i < (1 — a) for all / > 1 with /^ / = 1. By 
considering separately the positive and negative parts of a general / G T^(0,1) we 
deduce that ||T — 1 (g) < (1 — cl)- By dehnition Tq = q and T*1 = 1 so that 

(1 ® q)T = T{1 ®q) = l®q. This yields that {T — 1 ^ q)’^ = T^ — 1 ^ q. We thus 
obtain the hrst statement by noting that 

||f - 1 ® = ||(f - 1 <8) < (1 - a)K 

By taking adjoints in this inequality we see also that 

< (l-a)^ 


which clearly implies the second part of the Corollary. 


□ 
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The integrate-and-fire neuron 


We next apply the above notions to find unique periodic solutions to the periodic 
differential equations in our model and prove exponential convergence of any solution 
to the periodic one. Consider the linear and time-dependent differential equation in 


x'{t) = Q{t)x{t), (24) 

where x : [0, cxo) —)■ M” and Q : [0, oo) is a continuou^ matrix valued 

function that is 1-periodic. We assume that for each t the matrix Q is probability 
generating, i.e. its column sums are zero, diagonal elements are non-positive and 
other elements non-negative. 


Theorem 6.5. Assume that for some 0 < t^ < ti < 1 and for some constant 
c > 0 the 1-periodic probability generating matrix Q(t) = {qij) in (1241) satisfies the 
condition 

>c for t e {to, ti), i = 1 , 2 ,..., n, 

where we interpret gi,o =: qi,n- Then the system [2^ admits a unique 1-periodic 
solution (that is a probability distribution) and all solutions (that are probability 
distributions) converge with exponential rate towards this solution as t ^ oo. 

Proof. Denote by ^{t,s) the transition (operator) matrix of the flow map, so that 
x{t) = ^{t,s)x{s) for all 0 < s < t regardless of the initial values. We then know 
a priori that the coefficients of <h cannot be negative, since solutions with initial 
condition inside the positive cone stay inside the positive cone. For the proof of the 
lemma it is enough to show that 


a := 5{^{ti,to)) < 1. (25) 

Namely, then fl20|) and the equality *h(l,0) = <h(l, ti, )<F(ti, to)‘h(to, 0) shows that 

(5(<F(1,0)) < a < 1, 

whence Theorem 16.21 verifies that there is a unique probability vector (probability 
distribution on {1,2,.. .n}) that solves <F(l,0)7r = vr. This implies that there is a 
unique 1-periodic solution to the flow. Finally, for any integer m > 1 the periodicity 
of A yields that <h(m, 0) = (*h(l, 0))"*, whence the second statement of Theorem 16.21 
shows that ||x(m) — vr||i —)■ 0 with exponential rate and this clearly yields the stated 
exponential convergence of the flow to the periodic solution. By Theorem 16.31 it is 
enough to check that the transition matrix <F(ti,to) satisfies ^{ti,to)jy > 0 for all 
1 ^ j, k < n. Equivalently, we need to show that 


Xj{ti) > 0 for all j G (1, 2,..., n}, (26) 

where x solves the flow x' = Q{t)x on t G [to, ^i] with initial condition Xj{to) = 1 for 
j = k, and Xj{to) = 0 for j ^ k. By the cyclic symmetry of the situation we may 
as well assume that k = 1. Then the first equation and the assumed properties of 

^One can easily weaken the assumption of continuity 
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Q{t) yield the differential inequality (the constants c,C>0 below are fixed for all 
t G [to, ti]) 

x[ > —Cxi with xi(to) = 1 

which shows that 

xi > bi for to < t < ti, (27) 

where 6i > 0 depends only on C and ti — to- Then the last equation implies 

x'n > CXi — CXn > cbi — CXn with Xn(to) = 0 
and standard integration gives 

Xn(t) > bn(t - to) for to < ^ < h- 

Iterating, the differential inequalities x'f. > cxk-i — Cxk finally produce the bound 

Xk{t) > bk{t — to)"^^”^ for to < t < ti and 2 < k < n. (28) 

The statement fl26|l follows from (HT} and fl28|l . where one notes that the obtained 
bound depends only on n, c, C and ti — to- □ 

In the applications, the infinitesimal generator of the process is given by the 
matrix 


-s{t) s+{t) 



s~{t) —s{t) s~^{t) 



s-{t) 

-s{t) 

s^t) 

s^t) 

s-{t) 

-s+(t) 


with periodic and piecewise continuous functions s"*", s and s = s’*" + s . The 
transition mapping <h = $(0,1) then satisfies the required condition 

5(<h) < 1 

Indeed, the conditions of the previous theorem are clearly satisfied in our setup. 

Expectation values 

For a periodic density s(t), the expected time span between two consecutive firings, 
given that the neuron fired at Tq = a: is 

The expectation E{x) is a periodic function 

E{x + 1) = J (t — X — l)p{x + l,t)dt 
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The expected fc:th firing delay is calculated recursively from 


E,{x) 

Ek{x) 


E{x) 

x)pk{x,t)dt 

p{x, z)(t — z)pk-i{z, t)dzdt + 

j Ek-i{z)p{x,z)dz + Ei{x) 
T*Ek-i{x) + Ei{x) 



As a consequence 



x)p{x, z)pk-i{z, t)dzdt 


Ek{x) - Ek-i{x) = T*{Ek-i{x) - Ek-2{x)) 

= T*^-\E2ix) - E,{x)) 

= T*’‘-^E{x) 

As before, we shall denote by 'g the restriction of a 1-periodic function on M to the 
interval (0,1). By the following limit exists in L°°{0, 1) 

m{x) = lim {T*)^E{x) = lim {Ek{x) — Ek-i{x)), 

k^oo k^oo 


and the rate of convergence is exponential. Then rh{x) is a hxed point of T*, and 
Corollary 16.41 verihes that it is a constant function. Function m{x) is periodic hence 
it is a constant. 

Note that 

lim = lim 4t*^~^Ei{x) + ... + Ei{x)) 

fc—^oo k fc—>-00 k 

= m 


Both the differences Ek{x) — Ek-i{x) and the mean thus converge uniformly 
to the constant m. If g is a 1-periodic hxed point of T and hence of T, then 


(g,m) 



lim (g, f*^E) 

fc—>-oo 


lim (T*^g, E) 

k—yoo 


{g,E) 


The constant is then given as 

q(x)E(x)dx 

m = ——j- 

fo q{x)dx 

We have thus shown: 
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Theorem 6.6. Assume that E G L°°. Both the differences Ex — Ek- 
^ converge uniformly and exponentially towards the constant m. 

Remark. 

Theorem 14.21 shows that for the refractory neuron 


/ q{x)E{x)dx = 1 

Jo 

In the absence of a refractory period this is easily proved: 


E{x) = / {t — x)p{x,t)dt 

J X 

poo J 

Jx dt 

poo 


7 poo 

—E{x) = / s{x)e-J-^^^'^^^dt 

dx 


Since E{x) is periodic 


= —1 + s{x)E{x) 


d 

s(x)E(x)dx = dx+ —Eix)dx = 1 

./n ./n dx 


It has already been shown that q(t) = s(t) in this special situation. 


and the mean 
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